Interpretation for Iq-tree concordance factors at http://www.robertlanfear.com/blog/files/archive-2018.html

initialSet

First phylogenetic tree estimated with data from ortholog11-noBouvieri-queenconsambig.tar.gz, consensus alignments filtered through EAPhy, RAxML trees called as -s fasta -m GTRGAMMA -N 10 -p 12345 per locus and ASTRAL tree estimated (through script runRaxMLgeneTrees_ASTRALspeciesTree.R)

Showing tree comparison with published in Romiguier et al. 2017

Most samples

allMessor

mitogenome_iqtree

tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allMessor_busco_iqtree/concord.cf.tree")

tree_root <- root(tree, outgroup = "SRR4292898", edgelabel = TRUE)  ## by aciculatus as outgroup

nohyb <- dt %>%
    filter(label %in% tree$tip.label)

grp <- list()
for (sp in unique(nohyb$species)) {
    grp[[sp]] <- nohyb[nohyb$species %in% sp, "label"]
}

subgrp <- list()
for (sp in unique(nohyb$mitochondria)) {
    subgrp[[sp]] <- MRCA(ggtree(tree_root), nohyb[nohyb$mitochondrial %in% sp, "label"])
}
subgrp <- subgrp[names(subgrp)[!names(subgrp) %in% ""]]

grouped <- groupOTU(tree_root, .node = grp)
p <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    nohyb[, c("label", "caste")] + geom_tree(aes(color = group)) + geom_tiplab(size = 2,
    aes(color = group), hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
        breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
    "#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 28) +
    geom_nodelab(size = 1.8, nudge_x = 0.18)

for (nd in 1:length(subgrp)) {
    p <- p + geom_cladelabel(node = subgrp[[nd]], label = names(subgrp)[nd], offset = 2.2,
        align = TRUE, fontsize = 2)
}
grid.arrange(p)

busco_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allMessor_busco_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allMessor_busco_iqtree/geneTrees.treefile"))

pops3 <- left_join(data.frame(label = iqtree@phylo$tip.label), dt) %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(label, "_", lineage))

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w"), utf8ToInt("m")),
        breaks = c("queen", "worker", "male")) + scale_color_manual(values = c("#FFAD56",
    "#C5EABD", "purple"), breaks = c("queen", "worker", "male")) + xlim(0, 25) +
    geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

Using 2709 loci.

busco_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allMessor_busco_concatenatedRAxML/RAxML_bipartitionsBranchLabels.allMessor_busco')

tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
#                       .node= grp)

grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  ggnewscale::new_scale_color() +
  geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
                size = 2, na.rm = TRUE) +
  scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
                     breaks = c("queen", "male", "worker")) +
  scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
                     breaks = c("queen", "male", "worker")) +
  xlim(0, max(p$data$x)) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

allNoHybrid

mitogenome_concatenatedRAxML

tree <- read.raxml("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allNoHybrid_mitogenome_concatenatedRAxML/concatenated_raxml/RAxML_bipartitionsBranchLabels.mitogenome_75pConcatenated.tree")

dtTree <- dt %>%
    filter(label %in% tree@phylo$tip.label) %>%
    mutate(label2 = paste0(substring(hybrid, 1, 1), "_", label, "_", mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)) {
    grp[[sp]] <- dtTree[dtTree$species %in% sp, "label"]
}

# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in%
# 'SRR4292897')), .node= grp)

grouped <- groupOTU(tree, .node = grp)

p <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    dtTree[, c("label", "caste", "mitochondrial", "label2", "species")]

q1 <- p + geom_tree(aes(color = group)) + scale_color_manual(name = "species", values = spCol[names(spCol) %in%
    names(grp)]) + ggnewscale::new_scale_color() + geom_tiplab(size = 2, aes(label = label2,
    color = mitochondrial), hjust = -0.09, aligned = TRUE, show.legend = F) + scale_color_manual(values = linCol) +
    guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() + geom_tippoint(aes(x = x +
    0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) + scale_shape_manual(values = c(utf8ToInt("q"),
    utf8ToInt("m"), utf8ToInt("w")), breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
    "#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, max(p$data$x) +
    6) + geom_label2(aes(label = bootstrap, subset = as.numeric(sub("/.*", "", bootstrap)) >
    50), size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) + theme(legend.position = "bottom")

q1

busco_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allNoHybrid_busco_concatenatedRAxML/RAxML_bipartitionsBranchLabels.allNoHybrid_busco')

tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
#                       .node= grp)

grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  ggnewscale::new_scale_color() +
  geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
                size = 2, na.rm = TRUE) +
  scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
                     breaks = c("queen", "male", "worker")) +
  scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
                     breaks = c("queen", "male", "worker")) +
  xlim(0, max(p$data$x)) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

phasedH0_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allNoHybrid_phasedH0_iqtree/iqtreeOutput_Filtered_concord/phasedH0_min53IndFiltered_concord.cf.tree")

pops3 <- dt %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label = paste0(sra_access_number, ".0"), label2 = paste0(sra_access_number,
            "_", lineage), label3 = sra_access_number) %>%
    filter(label %in% iqtree@phylo$tip.label)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897.0")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
        breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
    "#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 21) +
    geom_nodelab(size = 2, nudge_x = 0.42)

allNoMaleHybrid

phasedH0_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allNoMaleHybrid_phasedH0_concatenatedRAxML/RAxML_bipartitionsBranchLabels.allNoMaleHybrid_phasedH0.tree')

tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
#                       .node= grp)

grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  ggnewscale::new_scale_color() +
  geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
                size = 2, na.rm = TRUE) +
  scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
                     breaks = c("queen", "male", "worker")) +
  scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
                     breaks = c("queen", "male", "worker")) +
  xlim(0, max(p$data$x)) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

AllNoMales

mitogenome_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/AllNoMales_mitogenome_concatenatedRAxML/concatenated_raxml/RAxML_bipartitionsBranchLabels.mitogenome_75pConcatenated.tree')

tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
#                       .node= grp)

grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  ggnewscale::new_scale_color() +
  geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
                size = 2, na.rm = TRUE) +
  scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
                     breaks = c("queen", "male", "worker")) +
  scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
                     breaks = c("queen", "male", "worker")) +
  xlim(0, max(p$data$x)) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

pureNoMales

mitogenome_iqtree

tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/pureNoMales_mitogenome_iqtree/iqtreeOutput_concord/pureNoMales_concord.cf.tree")

tree_root <- root(tree, outgroup = "SRR4292898", edgelabel = TRUE)  ## by aciculatus as outgroup

nohyb <- dt %>%
    filter(sra_access_number %in% tree$tip.label)

grp <- list()
for (sp in unique(nohyb$species)) {
    grp[[sp]] <- nohyb[nohyb$species %in% sp, "label"]
}

subgrp <- list()
for (sp in unique(nohyb$mitochondria)) {
    subgrp[[sp]] <- MRCA(ggtree(tree_root), nohyb[nohyb$mitochondrial %in% sp, "label"])
}
subgrp <- subgrp[names(subgrp)[!names(subgrp) %in% ""]]

grouped <- groupOTU(tree_root, .node = grp)
p <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    nohyb[, c("label", "caste")] + geom_tree(aes(color = group)) + geom_tiplab(size = 2,
    aes(color = group), hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
        breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
    "#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 20) +
    geom_nodelab(size = 1.8, nudge_x = 0.18)

for (nd in 1:length(subgrp)) {
    p <- p + geom_cladelabel(node = subgrp[[nd]], label = names(subgrp)[nd], offset = 2.2,
        align = TRUE, fontsize = 2)
}
grid.arrange(p)

messorQueens

mitogenome_iqtree

tree <- read.iqtree('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_mitogenome_iqtree/concord.cf.tree')

nLoci <- nrow(read.delim('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_iqtree/geneTrees.treefile'))

# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))

dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
                      .node= grp)
# tree2 <- tree
# tree2@phylo <- reorder(phangorn::midpoint(tree2@phylo)) ## to use a midpoint root

# grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
              hjust = -0.09,
              # align = TRUE,
              linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  ggnewscale::new_scale_color() +
  xlim(0, max(p$data$x)+0.1) +
  geom_treescale(x=0, y=0) + 
  geom_nodelab(aes(subset = as.numeric(sub("/.*", "", label)) > 50), 
               size = 2, hjust = 1.3, nudge_y = 0.15) +
  theme(legend.position = "bottom")

#   geom_label2(aes(label= label),subset = as.numeric(sub("/.*", "", label)) > 50)) + , size = 3, nudge_x = -0.5, label.padding = unit(0.1, "lines")
q1

mitogenome_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_mitogenome_concatenatedRAxML/RAxML_bipartitionsBranchLabels.messorQueens_mitogenome')

# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
                      .node= grp)

# grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  xlim(0, max(p$data$x)+0.1) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

busco_iqtree

tree <- read.iqtree('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_busco_iqtree/concord.cf.tree')

nLoci <- nrow(read.delim('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_iqtree/geneTrees.treefile'))

#tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))

dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
                      .node= grp)

# grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
              hjust = -0.09,
              # align = TRUE,
              linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  xlim(0, max(p$data$x)) +
  geom_treescale(x=0, y=0) + 
  geom_nodelab(aes(subset = as.numeric(sub("/.*", "", label)) > 50), size = 2, 
               hjust = 1.3, nudge_y = 0.15) +
  theme(legend.position = "bottom")

#   geom_label2(aes(label= label),subset = as.numeric(sub("/.*", "", label)) > 50)) + , size = 3, nudge_x = -0.5, label.padding = unit(0.1, "lines")
q1

busco_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_busco_concatenatedRAxML/RAxML_bipartitionsBranchLabels.messorQueens_busco')

# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
                      .node= grp)

# grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  xlim(0, max(p$data$x)) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

phasedH0_iqtree

tree <- read.iqtree('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_phasedH0_iqtree/concord.cf.tree')

# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))

dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
                      .node= grp)

# grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
              hjust = -0.09,
              # align = TRUE,
              linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
   xlim(0, max(p$data$x)+0.004) +
  geom_treescale(x=0, y=0) + 
  geom_nodelab(aes(subset = as.numeric(sub("/.*", "", label)) > 50), size = 2, hjust = 1.3, nudge_y = 0.15) +
  theme(legend.position = "bottom")

#   geom_label2(aes(label= label),subset = as.numeric(sub("/.*", "", label)) > 50)) + , size = 3, nudge_x = -0.5, label.padding = unit(0.1, "lines")
q1

phasedH0_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_phasedH0_concatenatedRAxML/RAxML_bipartitionsBranchLabels.messorQueens_phasedH0')

# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))

dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
                      .node= grp)

# grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
              hjust = -0.09,
              # align = TRUE,
              linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  xlim(0, max(p$data$x)+0.01) +
  geom_treescale(x=0, y=0) + 
  geom_nodelab(aes(subset = as.numeric(sub("/.*", "", label)) > 50), hjust = 1.3, nudge_y = 0.15, size = 2) +
  theme(legend.position = "bottom")

#   geom_label2(aes(label= label),subset = as.numeric(sub("/.*", "", label)) > 50)) + , size = 3, nudge_x = -0.5, label.padding = unit(0.1, "lines")
q1

Queen lineages subset

qsl

mitogenome_iqtree

tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qsl_mitogenome_iqtree/iqtreeOutput_concord/qsl_concord.cf.tree")

tree_root <- midpoint(tree, node.labels = "support")  ###midpoint root because for some reason forgot to include outgroup sample

nohyb <- dt %>%
    filter(sra_access_number %in% tree$tip.label)

grp <- list()
for (sp in unique(nohyb$species)) {
    grp[[sp]] <- nohyb[nohyb$species %in% sp, "label"]
}

subgrp <- list()
for (sp in unique(nohyb$mitochondria)) {
    subgrp[[sp]] <- MRCA(ggtree(tree_root), nohyb[nohyb$mitochondrial %in% sp, "label"])
}
subgrp <- subgrp[names(subgrp)[!names(subgrp) %in% ""]]

grouped <- groupOTU(tree_root, .node = grp)
p <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) + geom_tree(aes(color = group)) +
    geom_tiplab(size = 3, aes(color = group), hjust = 0, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp)], name = "species") + xlim(0, 12) + geom_nodelab(size = 1.8, nudge_x = 0.5)


for (nd in 1:length(subgrp)) {
    p <- p + geom_cladelabel(node = subgrp[[nd]], label = names(subgrp)[nd], offset = 1.7,
        align = TRUE, fontsize = 3)
}
grid.arrange(p)

phased_AstralMrBayes

path = "/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qsl_phased_AstralMrBayes/"
pops <- read.delim(paste0(path, "qsl_ID_pop.txt")) %>%
    rename(label = ID) %>%
    left_join(dt[, c("label", "species")])

tree <- read.tree(paste0(path, "/speciesTree/qslMrBayes_SpeciesTree.tree"))
grp <- list()
for (sp in unique(pops$species)) {
    grp[[sp]] <- pops[pops$species %in% sp, "label"]
}

pops2 <- pops %>%
    mutate(label2 = paste0(label, "_", lineage))

grouped <- groupOTU(root(tree, "SRR4292897"), .node = grp)
tree <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops2
tree$data$local.posterior <- ifelse(str_detect(as.numeric(tree$data$label), "[^a-zA-Z]"),
    tree$data$label, "")
tree1 <- read.tree(paste0(path, "/speciesTree/qslMrBayes_SpeciesTree_t1.tree"))
grouped1 <- groupOTU(root(tree1, "SRR4292897"), .node = grp)
tree1 <- ggtree(grouped1, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops2
tree$data$quartet.support <- ifelse(str_detect(as.numeric(tree1$data$label), "[^a-zA-Z]"),
    tree1$data$label, "")

tree8 <- read.astral(paste0(path, "/speciesTree/qslMrBayes_SpeciesTree_t8.tree"))
dt2 <- tree8@data %>%
    select(q1:q3, node) %>%
    drop_na()

pies <- nodepie(data.frame(dt2), cols = 1:3, color = c("red", "blue", "yellow"))
df <- tibble::tibble(node = as.numeric(dt2$node), pies = pies)
q <- tree %<+% df

library(ggpp)
q + geom_tiplab(aes(label = label2, color = group), size = 3, hjust = 0, aligned = TRUE) +
    geom_plot(data = td_filter(!isTip), mapping = aes(x = x, y = y, label = pies),
        vp.width = 0.1, vp.height = 0.07, hjust = 0.5, vjust = 0.5) + xlim(0, 12) +
    scale_color_manual(values = spCol[names(spCol) %in% names(grp)]) + geom_text2(aes(label = as.numeric(local.posterior),
    x = branch), size = 3, vjust = -0.2, color = "darkgreen") + geom_text2(aes(label = as.numeric(quartet.support),
    x = branch), size = 3, vjust = 1.1, color = "darkred")

phased_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qsl_phased_iqtree/phased_qsl.contree.tree")

iqtree <- treeio::drop.tip(iqtree, iqtree@phylo$tip.label[str_detect(iqtree@phylo$tip.label,
    "\\.1")])

pops3 <- pops %>%
    mutate(label2 = paste0(label, "_", lineage), label3 = label, label = paste0(label,
        ".0"))

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897.0")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3

treeIQ + geom_tree(aes(color = group)) + geom_tiplab(aes(label = label2, color = group),
    size = 3, hjust = 0, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + xlim(0, 13) + geom_nodelab(size = 3, nudge_x = 0.2)

- compare ASTRAL vs IQ-tree

T1 <- tree + geom_tree(aes(color = group)) + geom_tiplab(aes(label = label2, color = group),
    size = 2.5, hjust = 0, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + xlim(0, 15) + labs(title = "ASTRAL") + geom_nodelab(size = 3,
    nudge_x = 0.12) + theme(legend.position = "none", plot.margin = unit(c(0, 0,
    0, 0), "cm"))

d <- data.frame(old = pops3$label, new = pops3$label3)
groupedIQ@phylo = rename_taxa(groupedIQ@phylo, d, old, new)

pops32 <- pops3[, c("label3", "species", "label2")] %>%
    rename(label = label3)

treeIQ2 <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops32

T2 <- treeIQ2 + geom_tree(aes(color = group)) + geom_tiplab(aes(label = label2, color = group),
    size = 2.5, hjust = 1, aligned = TRUE) + scale_x_reverse(limits = c(13, 0)) +
    scale_color_manual(values = spCol[names(spCol) %in% names(grp3)]) + labs(title = "IQ-tree") +
    geom_nodelab(size = 3, nudge_x = 0.12) + theme(legend.position = "none", plot.margin = unit(c(0,
    0, 0, 0), "cm"))

d1 = T1$data[T1$data$isTip, c(1:2, 4:10)]
d1$x[] = 1
d2 = T2$data[T2$data$isTip, c(1:4, 7:11)]
d2$x[] = 2

TTcon <- rbind(d1, d2)

L1 = ggplot(TTcon, aes(x = x, y = y, colour = group, group = label)) + geom_line() +
    theme_void() + theme(legend.position = "none", plot.margin = unit(c(0, 0, 0,
    0), "cm")) + scale_color_manual(values = spCol[names(spCol) %in% names(grp3)])

cowplot::plot_grid(T1, L1, T2, nrow = 1, align = "hv", rel_widths = c(1, 0.5, 1))

qslAssemblies

qsl set of samples plus the three genome assemblies

busco_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_busco_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_busco_iqtree/geneTrees.treefile"))

pops3 <- left_join(data.frame(label = iqtree@phylo$tip.label), dt) %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(label, "_", lineage))

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
        "worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
    "worker")) + xlim(0, 10) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

Using 2559 loci.

phasedH0_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_phasedH0_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_busco_iqtree/geneTrees.treefile"))

pops3 <- left_join(data.frame(label = iqtree@phylo$tip.label), dt) %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(label, "_", lineage))

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
        "worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
    "worker")) + xlim(0, 10) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

Using 2559 loci.

phasedH01_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_phasedH01_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_phasedH01_iqtree/geneTrees.treefile"))

pops3 <- data.frame(ID = iqtree@phylo$tip.label) %>%
    mutate(label = word(ID, sep = fixed("."), 1, 1)) %>%
    left_join(dt) %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(ID, "_", lineage), label = ID)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
# groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in%
# 'SRR4292897')), .node= grp3)

groupedIQ <- groupOTU(iqtree, .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
        "worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
    "worker")) + xlim(0, 18) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

Using 3543 loci.

qslWorkerLineage

qsl set of samples + a worker per lineage

busco_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslWorkerLineage_busco_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslWorkerLineage_busco_iqtree/geneTrees.treefile"))

pops3 <- data.frame(ID = iqtree@phylo$tip.label) %>%
    mutate(label = word(ID, sep = fixed("."), 1, 1)) %>%
    left_join(dt) %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(ID, "_", lineage), label = ID)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
# groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in%
# 'SRR4292897')), .node= grp3)

groupedIQ <- groupOTU(iqtree, .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
        "worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
    "worker")) + xlim(0, 18) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

Using 2640 loci.

Queen lineages subset 2

Different samples chosen using only RNAseq dataset

sub1sp

Only one individual per species

VCF_SNAPP

beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2sub1sp_VCF_SNAPP/qsl2_vcfbedFiltered.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)
mean stderr_mean stdev variance median mode geom_mean hpd_interval_low hpd_interval_high act ess
posterior -10446.020 0.352 3.562 12.687 -10445.608 n/a n/a -10453.316 -10439.850 17609.569 102.274
likelihood -10418.099 0.271 2.889 8.344 -10417.808 n/a n/a -10424.397 -10413.449 15900.548 113.267
prior -27.921 0.191 2.032 4.127 -27.543 n/a n/a -32.237 -24.675 15967.076 112.795
lambda 0.179 0.003 0.058 0.003 0.172 n/a 0.169961734922052 0.077 0.292 3782.016 476.201
treeHeightLogger 10.985 0.186 1.765 3.114 10.825 n/a 10.8479409107085 8.050 14.773 20065.133 89.758
clockRate 0.002 0.000 0.000 0.000 0.002 n/a 0.00176233759024145 0.001 0.002 22665.337 79.461 
beast_log %>%
    pivot_longer(cols = -Sample) %>%
    ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()

beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2sub1sp_VCF_SNAPP/qsl2_vcfbedFilteredMCC.tree")

p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
    size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
    2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
    breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
    40, 6, 6))

qsl2

VCF_SNAPPsp

beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsp/qsl2_sp.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)
mean stderr_mean stdev variance median mode geom_mean hpd_interval_low hpd_interval_high act ess
posterior -7914.304 0.245 3.383 11.447 -7914.002 n/a n/a -7920.712 -7907.756 18896.626 190.563
likelihood -7884.929 0.238 2.993 8.957 -7884.626 n/a n/a -7890.772 -7879.766 22826.781 157.753
prior -29.375 0.066 1.516 2.297 -29.067 -28.4035981230664 n/a -32.592 -27.080 6883.582 523.129
lambda 0.272 0.002 0.073 0.005 0.266 0.276314482283504 0.262524665428786 0.145 0.424 2064.477 1744.267
treeHeightLogger 7.108 0.046 0.700 0.490 7.068 n/a 7.07354539315322 5.789 8.527 15243.478 236.232
clockRate 0.004 0.000 0.000 0.000 0.004 n/a 0.00372456049255618 0.003 0.005 13701.272 262.822 
beast_log %>%
    pivot_longer(cols = -Sample) %>%
    ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()

beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsp/qsl2_sp_MCC.tree")

p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
    size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
    2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
    breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
    40, 6, 6))

VCF_spSNAPP_AciMono

beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsp-AciMono/run1/qsl2_sp_AciMono.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)
mean stderr_mean stdev variance median mode geom_mean hpd_interval_low hpd_interval_high act ess
posterior -16831.287 0.226 3.447 11.880 -16830.815 n/a n/a -16838.153 -16825.185 3422.559 231.990
likelihood -16802.222 0.209 3.016 9.099 -16801.902 n/a n/a -16807.818 -16796.998 3801.807 208.848
prior -29.064 0.102 1.723 2.970 -28.671 n/a n/a -32.521 -26.495 2795.316 284.047
lambda 0.268 0.003 0.077 0.006 0.259 n/a 0.257946887799505 0.126 0.413 1460.433 543.674
treeHeightLogger 6.160 0.051 0.702 0.493 6.129 n/a 6.11966155006747 4.991 7.771 4110.355 193.171
clockRate 0.003 0.000 0.000 0.000 0.003 n/a 0.00332157547077767 0.003 0.004 4010.546 197.978 
beast_log %>%
    pivot_longer(cols = -Sample) %>%
    ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()

beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsp-AciMono/qsl2_sp_AciMonoCombinedMCC.tree")

p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
    size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
    2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
    breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
    40, 6, 6))

VCF_spSNAPP_noAci

beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-noAci_VCF_SNAPPsp/qsl2_sp_noAci.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)
mean stderr_mean stdev variance median mode geom_mean hpd_interval_low hpd_interval_high act ess
posterior -8099.867 0.139 3.117 9.718 -8099.446 n/a n/a -8106.267 -8094.486 8969.918 501.788
likelihood -8071.497 0.140 2.860 8.179 -8071.180 n/a n/a -8077.330 -8066.469 10722.871 419.757
prior -28.370 0.074 1.774 3.149 -28.123 n/a n/a -32.168 -25.422 7792.797 577.585
lambda 0.277 0.002 0.080 0.006 0.269 0.279568483126404 0.265861374445471 0.139 0.441 1676.720 2684.408
treeHeightLogger 6.980 0.034 0.684 0.467 6.931 n/a 6.94737794002518 5.626 8.284 11231.668 400.742
clockRate 0.003 0.000 0.000 0.000 0.003 n/a 0.00300838695349966 0.002 0.004 7267.469 619.335 
beast_log %>%
    pivot_longer(cols = -Sample) %>%
    ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()

beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-noAci_VCF_SNAPPsp/qsl2_sp_noAci_MCC.tree")

p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
    size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
    2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
    breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
    40, 6, 6))

VCF_spSNAPP_noAciAre

beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-noAciAre_VCF_spSNAPP/run1/qsl2_sp_noAciAre.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)
mean stderr_mean stdev variance median mode geom_mean hpd_interval_low hpd_interval_high act ess
posterior -16039.148 0.154 3.218 10.354 -16038.733 n/a n/a -16045.730 -16033.573 10264.242 438.513
likelihood -16013.579 0.135 2.805 7.866 -16013.271 n/a n/a -16019.087 -16008.492 10501.405 428.609
prior -25.569 0.074 1.562 2.440 -25.216 n/a n/a -28.756 -23.214 10214.810 440.635
lambda 0.322 0.002 0.090 0.008 0.313 n/a 0.310029883926007 0.165 0.506 1668.415 2697.771
treeHeightLogger 6.854 0.032 0.618 0.382 6.828 7.58032501811522 6.82618238315202 5.721 8.119 11876.655 378.979
clockRate 0.005 0.000 0.000 0.000 0.005 n/a 0.00475405329388956 0.004 0.006 12584.100 357.674 
beast_log %>%
    pivot_longer(cols = -Sample) %>%
    ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()

beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-noAciAre_VCF_spSNAPP/qsl2_sp_noAciAreCombinedMCC.tree")

p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
    size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
    2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
    breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
    40, 6, 6))

VCF_sampleSNAPP

beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsample/qsl2_sample.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)
mean stderr_mean stdev variance median mode geom_mean hpd_interval_low hpd_interval_high act ess
posterior -14875.395 0.480 5.584 31.181 -14874.934 n/a n/a -14886.017 -14864.935 18999.297 135.110
likelihood -14825.602 0.469 4.799 23.026 -14825.243 n/a n/a -14834.713 -14817.108 24508.675 104.738
prior -49.794 0.108 2.843 8.084 -49.388 n/a n/a -55.509 -45.099 3722.448 689.600
lambda 0.518 0.002 0.101 0.010 0.510 n/a 0.508478389893188 0.336 0.728 1355.567 1893.672
treeHeightLogger 8.200 0.026 0.766 0.587 8.163 n/a 8.16470357121235 6.667 9.600 3031.414 846.800
clockRate 0.003 0.000 0.000 0.000 0.003 n/a 0.0034195990990794 0.003 0.004 2850.903 900.416 
beast_log %>%
    pivot_longer(cols = -Sample) %>%
    ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()

beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsample/qsl2_sampleMCC.tree")

p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
    size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
    2), subset = posterior > 0.8), size = 4, nudge_y = 0.2) + scale_x_continuous(labels = abs,
    breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
    120, 6, 6))

mitogenome_iqtree

tree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_mitogenome_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_mitogenome_iqtree/geneTrees.treefile"))

pops3 <- dt %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
    filter(label %in% iqtree@phylo$tip.label)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
        "worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
    "worker")) + xlim(0, 15) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

busco_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_concatenatedRAxML/RAxML_bipartitionsBranchLabels.qsl2_busco_concatenated.tree')

tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
#                       .node= grp)

grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  ggnewscale::new_scale_color() +
  geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
                size = 2, na.rm = TRUE) +
  scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
                     breaks = c("queen", "male", "worker")) +
  scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
                     breaks = c("queen", "male", "worker")) +
  xlim(0, max(p$data$x)) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

mitogenome_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_mitogenome_concatenatedRAxML/RAxML_bipartitionsBranchLabels.qsl2_mitogenome')

tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
#                       .node= grp)

grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  ggnewscale::new_scale_color() +
  geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
                size = 2, na.rm = TRUE) +
  scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
                     breaks = c("queen", "male", "worker")) +
  scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
                     breaks = c("queen", "male", "worker")) +
  xlim(0, max(p$data$x)) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

busco_iqtree

tree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_iqtree/geneTrees.treefile"))

pops3 <- dt %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
    filter(label %in% iqtree@phylo$tip.label)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
        "worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
    "worker")) + xlim(0, 15) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

phasedH0_concatenatedRAxML

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_phasedH0_concatenatedRAxML/RAxML_bipartitionsBranchLabels.qsl2RNAseq_phasedH0')

tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
#                       .node= grp)

grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  ggnewscale::new_scale_color() +
  geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
                size = 2, na.rm = TRUE) +
  scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
                     breaks = c("queen", "male", "worker")) +
  scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
                     breaks = c("queen", "male", "worker")) +
  xlim(0, max(p$data$x)) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

phasedH0_concatenatedRAxML-2out

With two samples as outgroup

tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-2out_phasedH0_concatenatedRAxML/RAxML_bipartitionsBranchLabels.qsl2-2out_phasedH0')

tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
  filter(label %in% tree@phylo$tip.label) %>%
  mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))

grp <- list()
for (sp in unique(dtTree$species)){
  grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}

# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
#                       .node= grp)

grouped <- groupOTU(tree,.node= grp)

p <- ggtree(grouped, #branch.length='none', 
            ladderize=TRUE, right = TRUE)%<+%
  dtTree[,c('label','caste','mitochondrial','label2','species')] 

q1 <- p + 
  geom_tree(aes(color=group)) +
  scale_color_manual(name = "species",
                     values = spCol[names(spCol) %in% names(grp)]) +
  ggnewscale::new_scale_color() +
  geom_tiplab(size=2, aes(label = label2, color = mitochondrial), 
             hjust = -0.09,
             # align = TRUE,
                    linetype = "dotted", show.legend = F) +
  scale_color_manual(values = linCol) +
  guides(color=guide_legend(ncol=1)) +
  ggnewscale::new_scale_color() +
  geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
                size = 2, na.rm = TRUE) +
  scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
                     breaks = c("queen", "male", "worker")) +
  scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
                     breaks = c("queen", "male", "worker")) +
  xlim(0, max(p$data$x)) +
  geom_treescale() + 
  geom_label2(aes(label = bootstrap,
                  subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
              size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
  theme(legend.position = "bottom")

q1

Per clade

MbarDecCapOut

busco_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MbarDecCapOut_busco_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MbarDecCapOut_busco_iqtree/geneTrees.treefile"))

pops3 <- dt %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
    filter(label %in% iqtree@phylo$tip.label)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
        "worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
    "worker")) + xlim(0, 15) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

Using 2328 loci.

MbarQMWorkerLineage

busco_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MbarQMWorkerLineage_busco_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MbarQMWorkerLineage_busco_iqtree/geneTrees.treefile"))

pops3 <- dt %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
    filter(label %in% iqtree@phylo$tip.label)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
        breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
    "#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 11) +
    geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

Using 2263 loci.

MebeQ123Out

busco_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeQ123Out_busco_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeQ123Out_busco_iqtree/geneTrees.treefile"))

pops3 <- dt %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
    filter(label %in% iqtree@phylo$tip.label)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.035, shape = caste, color = caste), size = 3, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
        "worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
    "worker")) + xlim(0, 5) + geom_nodelab(size = 3, nudge_x = 0.3) + guides(color = guide_legend(ncol = 1))

Using 1654 loci.

MebeWasMinOut

busco_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeQ123Out_busco_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeQ123Out_busco_iqtree/geneTrees.treefile"))

pops3 <- dt %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
    filter(label %in% iqtree@phylo$tip.label)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
        "worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
    "worker")) + xlim(0, 5) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

Using 1654 loci.

MebeOut

amb_iqtree

iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeOut_amb_iqtree/concord.cf.tree")

nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeOut_amb_iqtree/geneTrees.treefile"))

pops3 <- dt %>%
    mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
        label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
    filter(label %in% iqtree@phylo$tip.label)

grp3 <- list()
for (sp in unique(pops3$species)) {
    grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
    .node = grp3)

treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
    pops3[, c("label", "label2", "species", "caste")]

treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
    size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
    names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
    geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
        breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
    "#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 13) +
    geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))

Using 1601 loci.

Ibericus

mitogenome_iqtree

tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/allIbericus_mitogenome_iqtree/concord.cf.tree")

nohyb <- dt %>%
    filter(sra_access_number %in% tree$tip.label)

p <- ggtree(tree, branch.length = "none", ladderize = TRUE, right = TRUE) %<+% nohyb[,
    c("label", "caste")] + geom_tree() + geom_tiplab(size = 2, hjust = -0.09, aligned = TRUE) +
    geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
        breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
    "#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 25) +
    geom_nodelab(size = 1.8, nudge_x = 0.4)

grid.arrange(p)

Ibericus queens

phasedH0_iqtree

tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/ibericusQueens_phasedH0_iqtree/concord.cf.tree")

nohyb <- dt %>%
    filter(sra_access_number %in% tree$tip.label)

p <- ggtree(tree, branch.length = "none", ladderize = TRUE, right = TRUE) %<+% nohyb[,
    c("label", "caste")] + geom_tree() + geom_tiplab(size = 3, hjust = -0.09, aligned = TRUE) +
    geom_tippoint(aes(x = x + 0.05, shape = caste, color = caste), size = 3, na.rm = TRUE) +
    scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
        breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
    "#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 10) +
    geom_nodelab(size = 1.8, nudge_x = 0.4)

grid.arrange(p)